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Abstract 

The prevalence of heterogeneous nucleation in nature was explained qualitatively by the classical 
theory for heterogeneous nucleation established over more than 60 years ago, but the quantitative 
validity and the key conclusions of the theory have remained unconfirmed. Employing the for¬ 
ward flux sampling method and the coarse-grained water model mW, we explicitly computed the 
heterogeneous ice nucleation rates in the supercooled water on a graphitic surface at various tem¬ 
peratures. The independently calculated ice nucleation rates were found to fit well according to the 
classical theory for heterogeneous nucleation. The fitting procedure further yields the estimate of 
the potency factor which measures the ratio of the heterogeneous nucleation barrier to the homo¬ 
geneous nucleation barrier. Remarkably, the estimated potency factor agrees quantitatively with 
the volumetric ratio of the critical nuclei between the heterogeneous and homogeneous nucleation. 
Our numerical study thus provides a strong support to the quantitative power of the theory, and 
allows understanding ice nucleation behaviors under the most relevant freezing conditions. 
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The freezing of water nearly all proceeds with the assistance from foreign substances, a 
process known as heterogeneous nucleation. In clouds, the dominant candidates for hetero¬ 


geneous ice nucleation are bacteria, pol 


molecular-weight organic compounds 


Q 


en grains, mineral dusts, soot particles, and high- 


2|. Despite its ubiquity, the microscopic picture 
behind such prevailing process still remains elusive, because of the complex and stochastic 
nature of the heterogeneous nucleation event. In particular, the mechanisms controlling 
heterogeneous ice nucleation are not well understood. 

Although a molecular understanding is still missing, the thermodynamic rationale behind 
the heterogeneous nucleation was already provided in 1950’s by the classical theory for het¬ 
erogeneous nucleation 3, 4|, an extension to the classical nucleation theory (CNT) for 
homogeneous nucleation, on the basis of macroscopic arguments. According to CNT, the 
formation of a critical nucleus needs to overcome a free energy barrier through spon¬ 

taneous fluctuations. In the case of a spherical solid nucleus forming from the supercooled 
liquid, the free energy barrier can be expressed as 
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where 7 /* is the solid-liquid interface free energy, Afiis is the chemical potential difference 
between liquid and solid, and p is the density of liquid. The homogeneous nucleation rate 
-Rhom varies with the nucleation temperature T following the Arrhenius equation jb]: 

Rhom = ^hom exp(-AGhoj^//cBT) , (2) 

where Ahom is the kinetic pre-factor. For homogeneous ice nucleation, both experiments [3- 
111 ] and simulations |l2l-ll7l| suggest that the temperature dependence of the homogeneous 
ice nucleation rate may be quantitatively described by CNT, with parametrization being 
rehned by the controlled experiments ISj. 

When a foreign flat wall (W) is present, the solid nucleus can preferentially form at the 
interface between the liquid and the wall (Fig. [11(a)). At its critical size, the solid embryo 
is under the unstable equilibrium with respect to the dissolution and the growth, which 
also indicates a mechanical equilibrium. At the liquid-solid-wall triple junction, solving the 
equation of equilibrium yields the Young’s equation 
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where and 'ysw are the surface tensions for liquid-wall and solid-wall interfaces. 6c defines 
the contact angle of solid embryo on the flat wall, with 6c = 0 and 6c = 180" indicating the 
complete wetting of the wall by solid and liquid, respectively. If the solid nucleus is further 
assumed to be part of the sphere, he., a spherical cap, its volume Kap can be expressed as 
heap /(^c) hsphere; where 

/(6'c) = (1 - COS 6'c)^(2 -h COS 6'c)/4 , (4) 


and ^sphere is the volume of the sphere containing the cap. Remarkably, under the framework 
of classical theory for heterogeneous nucleation, the factor f{6c) coincides with the ratio 
of the free energy barriers between the heterogeneous and homogeneous nucleation, he., 
f{6c) = AGhet/2^GLm- It thus follows that 
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Eqn. ([5]) provides a simple but robust explanation for the preference of the heterogeneous 
nucleation over the homogeneous nucleation: Instead of forming a spherical nucleus from 
spontaneous thermal fluctuations, only part of the sphere Rap needs to be nucleated when a 
foreign surface is present. Accordingly, the free energy barrier is reduced by the same factor 
f{6c) by which the volume of critical nucleus is reduced. Since /(R) measures the degree 
of the free energy reduction, it is also known as the potency factor. According to Eqn. (j4]), 
the potency factor f{6c) for a foreign wall is determined by the solid contact angle 6c, and 
varies between 0 and 1, as shown in Fig. [U^b). A wall with a lower solid contact angle yields 
a lower potency factor f{6c), thus further enhancing heterogeneous nucleation. Then the 
heterogeneous nucleation rate can be expressed by 


Rhet = Ahet exp{-f{6c)AGl^^/kBT) , 


( 6 ) 


Although the classical theory for heterogeneous nucleation offers a qualitative explana¬ 
tion to the prevalence of heterogeneous nucleation, its quantitative validity remains un- 
conhrmed. Auer and Frenkel 19| employed umbrella sampling method to calculate the 
nucleation barrier of the hard-sphere crystal that completely wets the smooth walls, and 
found that the computed barrier height is substantially higher than that predicted by the 
CNT. The disagreement was attributed to the omission in the CNT of the line tension at the 
liquid-solid-wall triple junction, which may become non-negligible when crystal completely 
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FIG. 1. (a) A solid critical nucleus forms at the interface between liquid and a foreign flat wall, 

with a solid contact angle Be- The contact angle Be can be determined by the surface tensions 
through Eqn. Q. The solid nucleus is assumed to be part of a sphere (dashed line) with a radius 
r*, i.e., a spherical cap with the volume Veap = /(^c)Kiphere- (b) The potency factor f{Bc) increases 
from 0 to 1 (homogeneous), as the solid contact angle Be varies from 0 to 180°. 


wets the wall. In a recent study by Winter et. al. 


20j, the total surface energies of the liquid 


nuclei (of both spherical and spherical cap shape) were directly obtained by Monte Carlo 
simulation for the Ising lattice gas model. It was found that the obtained surface energies 
could be comparable with the capillary approximation employed by CNT, if the line tension 
effects at the triple junction are considered. 


In this work, we show that the heterogeneous ice nucleation on a graphitic surface indeed 
supports the quantitative power of the theory. In particular, the validity of Eqn. ([5]) and (Ej) 
is strongly supported from the ice nucleation rates computed explicitly using the forward 
flux sampling method over a wide temperature range. Our work thus provides the hrst 
validation of heterogeneous CNT in the partially wetting regime where the potency factor 
is far enough from zero. 

Our molecular dynamics (MD) simulations were carried out using the coarse-grained 
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model of water (mW) 2l|. The inter-molecular interaction between water and carbon was 
adopted from a recent parameterization of the two-body term of the mW model, so that the 
strength of the water-carbon interaction reproduces the experimental contact angle (86°) 


of water on graphite 


22| . The model was recently employee 


study the heterogeneous ice nucleation on carbon surface 


22 


in direct MD simulations to 


24| . where the nonequilibrium 


freezing temperature of ice was found to increase due to the preferential nucleation of ice 


on carbon surface. Here we employ the forward flux sampling (FFS) method 


25 


26 | to 


systematically and explicitly compute the heterogeneous ice nucleation rates at various tem¬ 
peratures where spontaneous ice nucleation becomes too slow to occur in direct simulation. 
The details of the rate constant calculations can be found in the Supplementary Materials. 
Our MD simulation includes 4096 water molecules and 1008 carbon atoms, in a nearly cubic 
cell with a periodic boundary condition. The isobaric-isothermal canonical ensemble (NPT) 
with a Nose-Hoover thermostat was employed, with a relaxation time of 1 ps and 15 ps for 
temperature and pressure, respectively. A time step of 5 fs was used. It should be noted 
that while the homogeneous nucleation rate is measured by the nucleation frequency per 
unit volume, the heterogeneous nucleation rate should be characterized by the nucleation 
frequency per unit area. However because the simulation volume of liquid is small, and ice 
nucleation on carbon surface is strongly preferred, it is convenient to describe the heteroge¬ 
neous nucleation rate i?het on tho basis of volume, in order to facilitate a direct comparison 

with i?hom- 

Figure [2] shows the computed heterogeneous ice nucleation rates (in logarithm) as a 
function of nucleation temperature, in the range of 220 K to 240 K. To quantitatively explore 
the catalytic activity of the graphitic surface, we compare the obtained heterogeneous ice 
nucleation rates with the reported homogeneous ice nucleation rates from the previous work 
12| using the FFS method and the mW water model, as shown in Figure El It is clear 
that a graphitic surface yields the signihcantly enhanced ice nucleation rates, under all 
temperatures studied. Our results thus support the hnding by Lupi at. al. |22| and conhrm 
the enhanced ice nucleation capacity of carbon surface. 

The calculated heterogeneous ice nucleation rates at various temperatures allow assessing 
the quantitative validity of Eqn. ([6]). To do this, we £t the obtained heterogeneous ice 
nucleation rates at various temperatures according to the theory of nucleation, u sing the 


procedure employed by Li at. al. in analyzing the homogeneous ice nucleation 


12j. In 
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FIG. 2. Temperature dependence of ice nucleation rate (logarithm) in the mW water model, for 
both homogeneous ice nucleation (red) and the heterogeneous nucleation on a graphitic surface 
(blue). The solid red dots and blue squares represent the calculated ice nucleation rates by using 
the FFS method. The dashed red and blue lines indicate the fitting, on the basis of the CNT for 
homogeneous nucleation (Eqn. ([7|)) and the extension for heterogeneous nucleation (Eqn. ([8|)), 
respectively. The data of homogeneous nucleation were extracted from Ref. 


this procedure, the chemical potential difference Ap is approximated as a linear function of 
temperature, i.e., A/i = H{T — Tm)/Tm, where is the equilibrium melting temperature 
(274.6 K) of ice in the mW model, and is a constant; The liquid-solid interface energy 
'yis is assumed to be temperature independent. It is noted that both assumptions have been 
verihed for the mW water model in different simulation studies [^, 27|. For homogeneous 
ice nucleation, it was shown that the independently calculated homogeneous ice nucleation 
rate can be fitted according to the following expression: 


ln(/2hom) — lll(^hom) T 




hom 


(T - Tm^T ’ 


( 7 ) 


where In(Ahom) = 114.07±1.86 and Ghom = H‘^) = —3.72 ±0.08 x 10^ 


are the fitting constants 12|. We note that the nucleation barrier AG^om = ksChom/(T — 
TjnY- The htting yielded an estimate of = 31.01 ± 0.21 mJ m“^, which agrees well with 
the surface tensions computed through other approaches for the mW water model 17|, |28 |. 

Remarkably, the obtained heterogeneous ice nucleation rates are also found to ht well 
the classical theory for heterogeneous nucleation, as shown in Figure [2l Specifically, the 
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calculated ice nucleation rates -Rhet can be fitted according to 


ll^(-Rhet) ll^(^het) 


c, 


het 
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The fitting yields the estimate of the kinetic pre-factor for heterogeneous nucleation 
ln(y4het) = 102 ± 7.70, which is consistent with that for the homogeneous nucleation 
In(Ahom) = 114.07 ± 1.86 12|. More importantly, the other htting constant Chet = 
— 1.70 ± 0.07 X 10^ allows estimating the reduction of the nucleation barrier, as 
Chet/Chom = = /^(^c)- By Comparing the htting constants C from the het¬ 

erogeneous and homogeneous ice nucleation, we obtain the potency factor for the graphitic 
surface fbiOc) = 0.456 ± 0.019. This corresponds to a solid contact angle of 6c ~ 86 . 6 °. 
It is noted here that 6c is the contact angle between ice and graphene, and should not be 
confused with the water-graphene contact angle 86 °, although their magnitudes coincide 
here. 

It is then of interest to further test the validity of Eqn. ([5]), namely, the potency factor 
fb{6c) can be also quantitatively related to the volumetric ratio fv{6c) of the critical nu¬ 
cleus of the heterogeneous nucleation and the homogeneous nucleation. We note that such 
verihcation becomes possible in our study because the size of the critical nucleus A* can 
be independently estimated from the ensemble of nucleation trajectories obtained in the 
FFS calculation. Using the dehnition that the critical nucleus has the equal probabilities 
of dissolving and growing completely, i.e., with a committor probability pg = 0.5 j^, we 
obtained the estimate of the critical nucleus size A (number of water molecules contained in 
the critical ice nucleus) at various nucleation temperatures, as shown in Fig. |3l According 
to CNT, the critical size of the spherical nucleus in homogeneous nucleation is expressed by 
Ahom = 327r7^/(3p^(A/i)^). For mW water model, since 7 is nearly temperature indepen¬ 
dent and Ap = H{T — Tm)/Tm jl^, the critical size A^^^^ exhibits the following temperature 


dependence: 


KoraiT) = 


327r7^ 
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where i?hom is the temperature independent constant. The obtained critical size at 
various temperatures are found to £t well according to Fan fl9ll. as shown in Fig. [3](a). 


The good £t is not unexpected because previous studies [l3|, ll5| have shown the critical 


ice nucleus from homogeneous nucleation is nearly spherical. The fitting procedure yields 
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FIG. 3. (a) Variation of the critical ice nucleus size with temperature for both homogeneous 

and heterogeneous nucleation. The critical size is estimated based on the committor probability 
analysis, which yields 145 ± 5 at 225 K, 205 ± 10 at 230 K, 288 ± 7 at 235 K, and 410 ± 
10 at 240 K, respectively (See Supplementary Materials for more details). The critical size for 
homogeneous nucleation was obtained in Ref. [l^. The simulation data and the fitted curves are 
represented by data points and solid lines, respectively, (b) A snap shot of the critical ice nucleus 
forming from homogeneous ice nucleation at 240 K. (c) A snap shot of the critical ice nucleus 
forming on graphene surface at 240 K. 

the constant i?hom = 1-752 ± 0.027. For the heterogeneous ice nucleation, the same htting 
procedure was found to equivalently apply to the calculated critical nucleus size AJet, through 



( 10 ) 


which yields the constant i?het = 0.840 ± 0.014. By comparing the two htting constants 
fi, one obtains the volumetric ratio fv{0c) = AJet/Ahom = -Bhet/.Bhom- Remarkably, the 
obtained volumetric ratio fv{0c) = 0.480 ± 0.011 agrees quantitatively with the potency 
factor fb{0c) = 0.456 ± 0.019 estimated from the nucleation barriers. The quantitative 
validity of Eqn. ([5]), an important conclusion from CNT and its extension, is thus strongly 
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supported through the molecular simulation results based on the mW water model. 

The verihed quantitative validity of the CNT (and its extension) then allows predicting 
ice nucleation behavior in the presence of a heterogeneous nucleation center. The nucleation 
efficacy of the foreign surface can be generally described based on its potency factor /. Using 
the htted kinetic pre-factor Ahom (~ ^het) and Uhom, and Eqn. ([8]), one obtains both the 
heterogeneous ice nucleation rate and the corresponding critical nucleus size, as a function of 
the nucleation temperature, for different potency factors /. As shown in Fig. 01 the predicted 
ice nucleation rates clearly indicate the preference and the relevance of the heterogeneous 
nucleation at the moderate and low supercooling. For example, at 250 K, a nucleation 
center with the potency factor of / = 0.1 (equivalent to a solid contact angle 9c = 52.5°) 
yields an ice nucleation rate about 100 orders of magnitude higher than that of homogeneous 
nucleation at the same temperature. Intriguingly, the sizes of the critical nuclei from these 
relevant nucleation events fall within the range of a few hundred to a few thousand water 
molecules. This implies that the most relevant ice nucleation events mediated by an effective 
nucleation center under a low supercooling, i.e., —10°C ~ —20°C, can be possibly modeled 
by molecular simulations through using a reasonable number (10^ ~ 10^) of water molecules. 

Interestingly, when ice nucleation rate R is fixed, the theory of nucleation predicts that the 
critical size A* decreases with the potency factor / (Fig. 01(c)), and a homogeneous nucleation 
yields the minimum critical nucleus. This may appear surprising, but the prediction can be 
understood by the fact that the heterogeneous nucleation producing the same ice nucleation 
rate occurs at a much elevated temperature (Fig. 01(a)). For example, an ice nucleation 


rate of 10"^^ m“'^s“^ would require a homogeneous nucleation temperature Thom ~ 225 K, 
but a heterogeneous nucleation temperature Thet ~ 260 K for the nucleation center with 
a potency factor / = 0.1. The prediction (Fig. 01(c)) shows the critical nucleus for such 
heterogeneous nucleation at 260 K contains about 1220 water molecules {i.e., 1/10 of the 
critical size of the homogeneous nucleation at 260 K), larger than the critical size (310) 
of the homogeneous nucleation at 225 K. As the solid contact angle 9c decreases with the 
nucleation efficacy (Fig. [T]), a strong ice nucleation center yields a more “flat” ice nucleus 
that appears increasingly two-dimension like. It is noted that in such scenario the possible 
effect from the line tension at the triple junction can be non-negligible 19|, l20|. However 
its quantitative effect in ice nucleation rate is unclear. Using the density of ice, one can 
estimate the radius of the spherical segment {i.e., the frustum of the spherical cap) to be of 
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FIG. 4. The predicted temperature variation of (a) ice nucleation rate R{T) and (b) critical size 
of ice nucleus A*(T) in the mW water model, due to the presence of a heterogeneous nucleation 
center with different potency factor /. A small potency factor / indicates a strong ice nucleation 
efficacy, and / = 1 corresponds to homogeneous ice nucleation. (c) The predicted variation of the 
critical size A* for different potency factor /, subject to a fixed nucleation rate. 


the order of a few nano meters, implying that the dimension of an effective ice nucleation site 
is typically of the order of 10^ nm. This estimate may be used as an important parameter in 
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experiments for potentially observing ice nucleation in situ and designing effective strategy 
for controlling ice nucleation. 
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CALCULATION OF ICE NUCLEATION RATES 


Q 


We employed the forward flux sampling (FFS) method [ij to compute the ice nucleation 
rate on graphitic surface. The method has been applied to successfully study the homo¬ 


geneous ice nucleation [2[ using the mW water model. In FFS, the rate constant for the 
transition from the basin A to the basin B is obtained through the “effective positive flux” 
expression j^: Rab = ®Ao-P('^b|'^o)) where <Fao and P{Xb\Xo) ai’e the initial flux rate crossing 
the hrst interface Aq from basin A, and the probability for a trajectory starting from Aq to 
reach basin B, respectively. The interfaces are dehned based on an order parameter A, i.e., 
number of water molecules contained in the largest ice cluster. The ice-like and liquid-like 
water molecules are identihed based on the local order parameter gg Q , and an ice-like water 
molecule is dehned as the one with gg > 0.5. More details of computing gg can be found in 
Ref. { 2 ]. The initial hux rate can be computed by direct molecular dynamics simulation, 
through <Fao = No/itoV), where Nq is the number of successful crossings to the interface Aq 
from basin A, to is the total time of the initial sampling, and V is the simulation volume. 
As our recent work ^ has demonstrated the importance of the initial sampling on the relia¬ 
bility of the hnal nucleation rate, particularly in a heterogeneous environment, a sufficiently 
large sampling time {to > 300 ns) is used for computing $Ao to ensure the convergence. 
The growth probability P{\b\Xo) is obtained through P{Xb\\o) = 0^=1the 
individual crossing probability P(Aj|Ai_i) = where W is the number of successful 

crossings (~ 110) to the interface Aj, and Mj_i is the number of trial shootings (typically 
of 10^) from the interface Aj_i. The more details for computing P{Xb\Xo) are explained in 
Ref. |2|, [SJ . The statistical uncertainty of nucleation rate is mainly attributed to the error 
in the calculated P(Ai|Ai_i), which includes the variance of the binomial distribution of W 
and the landscape variance of the conhgurations collected at the previous interface Aj_i [l|. 
The calculated initial flux rate ^Ao, growth probability P(As|Ao), and the hnal nucleation 
rates R are listed in Table [I] for the heterogeneous ice nucleation at 225 K, 230 K, 235 K, 
and 240 K. The computed growth probability P(A|Ao), as a function of nucleus size A, is 
also shown in Fig. [1] 


At 220 K, ice nucleation becomes accessible in direct molecular dynamics simulation. 
Three out of four trajectories lead to spontaneous ice crystallization, with a total simulation 
time around 30 ns, which yields an estimate of the nucleation rate of the order of 10^^ m“^s“^. 
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FIG. 1. Calculated growth probability P(A|Ao) as a function of ice nucleus size A for heterogeneous 
ice nucleation on a graphitic surface at 225 K, 230 K, 235 K, and 240 K. The interfaces Aj used in 
FFS simulation correspond to the abscissa of the data points. 

TABLE I. Calculated initial flux rate growth probability P(Ab|Ao), and nucleation rate R 
for heterogeneous ice nucleation on a graphitic surface, using the mW water model and the FFS 
method. The critical nucleus size A* is determined by the committor probability analysis, on the 
basis of the computed growth probability. 


T/K 

^Ao/i 

m' 

-3g-l 

^AbIAo) 


Nucleation 

rate 

R 

/m 3s 1 

Critical size A' 

225 

1.69 

X 

1034 

7.64 ±2.16 X 

10’ 

-4 

1.29 

± 

0.36 

X 

1034 

145 ±5 

230 

1.12 

X 

1034 

4.83 ±1.83 X 

10- 

-6 

5.41 

± 

2.05 

X 

1028 

205 ± 10 

235 

8.90 

X 

1033 

1.51 ±0.74 X 

10- 

-10 

1.34 

± 

0.66 

X 

1024 

288 ±7 

240 

4.19 

X 

1033 

2.23 ±1.49 X 

10- 

-15 

9.34 

± 

6.24 

X 

1048 

410 ± 10 


The size of the critical ice nucleus A* can be numerically determined by the committor 
probability pg jo]. A critical crystalline nucleus is dehned as the one with 50% of probability 
of growing completely into crystal, i.e., pg = 0.5. The committor pg as a function of nucleus 
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Number of water molecules in ice nucleus 


FIG. 2. Calculated committor pg as a function of the nucleus size A for heterogeneous ice nucleation 
on a graphitic surface at various nucleation temperatures. The dotted lines indicate the upper and 
lower boundaries of the calculated pg, on the basis of the estimated error bar of pg. The horizontal 
dash line corresponds to a pg = 0.5, and intersects the both boundaries of the calculated pg, which 
yields the estimate of the uncertainty of A*. 


size A is obtained on the basis of the compnted growth probability P(Ai|Ai_i), and the critical 
size A* nnder different temperatnres are determined (Fig. [2]) and listed in Table [Tl 

To nnderstand whether the size of the nnclens A is a good approximation to the actnal 
reaction coordinate, we carried ont the committor distribntion analysis near the critical size 
A* ~ 205 for the heterogeneous nucleation at 230 K. 110 conhgurations were selected from the 
interface A = 210 and each one receives 20 shootings with Boltzmann distributed velocities. 
The obtained committor distribution is shown in Fi^El along with the Gaussian distribution 
with the intrinsic mean /i and standard deviation a [7| . The calculated distribution is peaked 
at /i = 0.521, with an overall agreement with the intrinsic committor distribution, indicating 
that A fairly accurately describes the actual reaction coordinate, and is thus a good order 
parameter. 
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FIG. 3. The computed histogram of committor pg for heterogeneous ice nucleation at 230 K at the 
interface A = 210. The Gaussian distribution function with the same intrinsic mean and standard 
deviation is shown as the red line. 
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